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Direct numerical simulation of turbulent 
flow over a backward-facing step 


By H. Le AND P. Moin 


1. Objectives 

The objectives of this study are: (a) to conduct a direct numerical simulation 
of turbulent backward facing step flow using inflow and outflow conditions and 
(b) to provide data in the form of Reynolds stress budgets for Reynolds averaged 
modeling. The report presents the basic statistical data and comparisons with the 
concurrent experiments of Jovic and Driver (1991, 1992) and budgets of turbulent 
kinetic energy. 

2. Accomplishments 

2 A Method 

The Navier-Stokes equations are discretized using a finite difference method on 
a staggered mesh. Uniform mesh spacing is applied in the streamwise (x) and 
spanwise (z) directions. In the vertical (t/) direction, non-uniform mesh is employed 
with mesh refinement at the wall and near the step. The fractional step method from 
Le and Moin (1990) is used for time advancement. The Navier-Stokes equations 
are first advanced using a second-order semi-implicit method without the pressure 
terms. The pressure is calculated by solving the Poisson equation, and the velocities 
are then corrected to satisfy the continuity equation. 

The mean velocity profile obtained from a boundary layer simulation by Spalart 
(1986) is imposed at the inlet at Rte = 667. Random velocity fluctuations u * , v\ 
and w r are superimposed on this profile according to a variant of the method of 
Lee et al (1992). The fluctuations are prescribed such that, at the inlet, the tur- 
bulence intensities and Reynolds shear stress of Spalart’s data are also duplicated. 
A convective boundary condition (Pauley et a/., 1988) is imposed at the exit. 

The streamwise domain consists of an entry section of length 10 h prior to the 
step and a 20 h post-expansion section, where h is the step height. The vertical 
dimensions before and after the expansion are W\ — 5 h and W 2 = 6 h which give an 
expansion ratio ER of 1.20. The spanwise dimension is 4 h where periodic bound- 
ary conditions are imposed. The simulation uses 770 x 194 x 66 grid points in the 
streamwise, wall normal, and spanwise directions, respectively. The Reynolds num- 
ber, based on h and the mean inlet free stream velocity is Rek = 5100. The 
computation uses 13 megawords of memory and requires approximately 55 CPU 
seconds per time step on a single processor CRAY Y-MP at a rate of 186 mflops. 
Statistical quantities axe averaged over the spanwise direction and time. About 
1100 CPU hours were required to obtain an adequate statistical sample. The total 
computational time corresponds to approximately 4.5 flow-through times. 
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FIGURE 1. Instantaneous spanwise vorticity contours, lj z \ min = —7, max = 4.95, 
increment = 1.3. 


2.2 Results 

Figure 1 shows contour plots of the instantaneous spanwise vorticity w* on a typi- 
cal vertical plane. The vorticity is normalized by Uo/h. A free-shear layer spreading 
from the step and interacting with the lower wall near the mean reattachment lo- 
cation, xr = 6 h, is discernible. 

The basic statistical quantities are compared to results from concurrent exper- 
iments by Jovic and Driver. In 1991, they conducted a backward facing step ex- 
periment at Rek = 6800 and ER = 1.09, herein referred as “JD1”. However, the 
results indicated that these parameters are not sufficiently close to those used in 
simulations. Thus, Jovic and Driver in 1992 conducted a second experiment with 
Reh = 4950 and ER = 1.20 (“JD2”). In the following sections, JD2 results will be 
used for comparison with numerical simulations. 

2.2.1 Reattachment length 

The instantaneous velocity fields indicate that the reattachment location oscillates 
in the streamwise direction and time and also varies in the spanwise direction. 
Several methods were used to determine the mean reattachment location, xr: (a) 
by the location at which the mean velocity U — 0 at the first grid point away from 
the wall, (b) by the location of zero wall-shear stress ( t w = 0), and (c) by the 
location of the mean dividing streamline (t/> = 0). A pdf method was also used in 
which the mean reattachment point is indicated by the location of 50% forward flow 
fraction. The pdf method was also used experimentally by Westphal et al. (1984) 
and Adams ef al. (1984). The results of the first three methods are within 0.1% of 
each other and about 4% from the pdf result. The calculated mean reattachment 
length is 6.0/( compared to 6.1/i measured in JD2 experiment. 

2.2.2 Pressure distribution 

The streamwise pressure coefficient at the step-wall is compared with the JD2 
results in Fig. 2. The pressure coefficient is normalized by the mean inlet free- 
stream velocity. The two symbols in Fig. 2 are for the two walls of the double-sided 
expansion in the JD2 experiment. The comparison shows an excellent agreement 
between computation and experiment. 
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FIGURE 2. Step- Wall pressure coefficient as a function of x. , computation; 

o , Jovic &; Driver (1992), bottom wall; a , Jovic & Driver (1992), top wall. 
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FIGURE 3. Step-wall coefficient of friction, 
(1992). 


, computation; • , Jovic &; Driver 
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FIGURE 4. Mean velocity profiles. , computation; • , Jovic &; Driver (1992). 

2.2.3 Skin friction coefficient 

The skin friction coefficient Cf is computed from the mean velocity profile and 
shown in Fig. 3. Also plotted are the Cf data from the JD2 experiment. The large 
peak negative skin friction in the recirculation region (~ —3 x 10“ 3 ), seen in both 
computation and experiment, is nearly 3 times larger than previous experimental 
findings, e.g. Adams et al (1984). Jovic and Driver (1991) also measured the skin 
friction in the reverse flow region for a wide range of Reynolds numbers. Their 
results show that the recirculation Cf decreases with increasing Reynolds number. 
The peak negative C / reaches a value of approximately — 1.0 x 10" 3 as Ren ~ 40000 
which is the Reynolds number range used in many experiments. Thus the large skin 
friction in the current study is due to low Reynolds number effects. 

2.2.4 Mean velocity profiles 

Figure 4 shows the mean velocity profiles compared with JD2 experimental results 
at three x-locations: x/h = 4 (recirculation region), x/h = 6 (reattachment loca- 
tion), and x/h = 10. Again, excellent agreement is obtained between computation 
and JD2 experiment. 

The near- wall mean velocity profile in the recovery region at x/h = 19 is compared 
to experimental data in Fig. 5. The shear velocity u T of the JD2 experiment is 
calculated from the skin friction coefficient of Fig. 3. Above = 10, both profiles 
show a noticeable shift below the log-law, indicating that the turbulent boundary 
layer has not fully recovered. Previous studies, however, reported a recovery of 
the log-law profile even as early as x — xr — 6 ft, e.g. Westphal et al. (1984). 
The apparent discrepancy is attributed to the method of obtaining the wall-shear 
velocity u T . In these experiments, the wall-shear velocity was calculated using the 
Clauser chart with the inherent assumption that the log-law of the zero-pressure 
gradient turbulent boundary layer is applicable. The result was lower values of u T . 
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Figure 5. Mean velocity profile at xfh = 19. , = y + ; , 

— (l/.41)log (?/ + ) + 5.0; , computation; • , Jovic & Driver (1992). 


Direct C / measurements using a pulsed wall probe by Westphal et al. confirm that 
the correct shear velocity in this flow region is indeed higher than that predicted 
by the Clauser chart. In the present configuration, this discrepancy is as high as 
17%, i.e., u T /u Tc « 1.17, where the subscript c denotes the value obtained using the 
Clauser chart. The deviation of the velocity profile from the log-law may be due to 
the effect of the strong streamwise adverse pressure gradient which is experienced 
by the flow following the sudden expansion (Nagano et a/., 1991). It may also be 
due to non-equilibrium effects which are persistent after reattachment. 

2.2.5 Turbulence intensities and Reynolds stress 

The turbulence intensity and Reynolds shear stress profiles, normalized by the 
inlet free-stream velocity Uo y are shown in Figs. 6-8. They are compared with 
those from JD2 experiment data at streamwise locations, x/h = 4, 6, 10. The 
overall agreement is good. 

2.2.6 Turbulent kinetic energy budget 

The terms of the turbulent kinetic energy transport equation are evaluated from 
the equation 


dt 
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FIGURE 6. Turbulence intensity normalized to inlet mean velocity Uq. 

, computation; • , Jovic &: Driver (1992), 




FIGURE 7. Turbulence intensity normalized to inlet mean velocity Uq. 

, computation; • , Jovic Sz Driver (1992). 
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-uV/6 r 0 2 

FIGURE 8. Reynolds shear stress — itV, normalized to inlet mean velocity Uq. 
, computation; • , Jovic & Driver (1992). 

The terms on the right hand side axe identified as follows: 

C k = Convection, 

Pk = Production, 

Tk — Turbulence Diffusion, 

D k = Viscous Diffusion, 
e k = Viscous Dissipation, 

II k = Velocity Pressure-Gradient. 

The budget for the turbulent kinetic energy, ^ is shown in Figs. 9 

and 10. 

At two step heights before the separation, the energy budget (Figs. 9(a) and 
10(a)) is similar to that of a turbulent boundary layer (Spalart, 1988) although 
there is an enhancement of the viscous terms near the wall. The budget for the 
recirculation region is shown in Fig. 9(b). This energy budget is very similar 
to that of a plane mixing layer (Bradshaw and Ferriss, 1965). This budget also 
agrees qualitatively with the measurement by Chandrsuda and Bradshaw (1981) 
for a backward facing step flow. Both production Pk and viscous dissipation e* 
have maximum values at the same point in the free- shear layer. The peak e k is 
approximately 60% of the production peak. The other 40% is balanced by the 
turbulence transport T*. Far from the wall, P* is the only major gain term. T k is 
a consuming term for 0.3 < yjh < 1, but a gain term outside of this range. Above 
the step {yjh > 1), T k is in balance with C*, the turbulence convection by the 
mean flow. As one approaches the wall, the production becomes a consuming term 
because of the negative gradient of the mean reverse flow (Fig. 10(b)) although its 
magnitude is relatively small. Very close to the wall, the two viscous terms, Dk and 
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FIGURE 9. Budget terms for turbulent kinetic energy, \ away from 

the wall, normalized by U^/h. (a) xjh = —2.0; (b) xjh = 4. , production; 

, turbulence transport; , viscous diffusion; , viscous dissipation; 

, velocity pressure-gradient; , convection. 










FIGURE 10. (cont.) Budget terms for turbulent kinetic energy near wall, normal 
ized by ui/v. (c) x/h = 10; (d) x/h = 18. For caption see previous page. 
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ejt, grow rapidly due the high intensities in all three directions (their value at the 
wall is about 40% higher than the peak production in the shear layer). The velocity 
pressure- gradient term becomes significant only at y + < 8 where it balances the 
viscous terms plus T* . The profiles have similar shapes through the reattachment 
point. 

All terms decay with x; however, the energy in the shear layer decays much faster 
than at the wall (the streamwise decay of the budget terms are not apparent in Fig. 
10 because they are normalized with the local u r ). By x/h = 10, the value of Dk at 
the wall is about 3 times the peak production value. Near the flow exit, y/h = 18, 
the turbulent kinetic energy budget resembles that of a boundary layer. However, 
the effects of the free-shear layer is still apparent, e.g., Tk is still large at yjh = 1 
and the peak tk is only ~ 85% of the peak Pk- 

3. Conclusions 

A direct numerical simulation of turbulent flow over a backward facing step at 
Re h — 5100 was successfully completed. There is good agreement in turbulence 
statistics deduced from simulations and concurrent experiments of Jovic and Driver. 
Of interest are two observations not previously reported for the backward-facing 
step flow: (a) at the relatively low Reynolds number considered, large negative skin 
friction is seen in the recirculation region; the peak |C/| is about three times the 
value reported at high Reynolds numbers; (b) the velocity profiles in the recovery 
region fall below the universal log-law. 

A large data base from this recently completed simulation has been archived. It 
contains up to third-order statistics at all locations in the recirculation, reattach- 
ment, and recovery zones and the budgets of all components of the Reynolds stress 
tensor. 
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